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ABSTRACT 


Observations revealed rich dynamics within prominences, the cool (10^ K), 
macroscopic (sizes of order 100 Mm) ‘clouds’ in the million degree solar corona. 
Even quiescent prominences are continuously perturbed by hot, rising bubbles. 
Since prominence matter is hundredfold denser than coronal plasma, this bub¬ 
bling is related to Rayleigh-Taylor instabilities. Here we report on true macro¬ 
scopic simulations well into this bubbling phase, adopting a magnetohydrody¬ 
namic description from chromospheric layers up to 30 Mm height. Our vir¬ 
tual prominences rapidly establish fully non-linear (magneto)convective motions 
where hot bubbles interplay with falling pillars, with dynamical details includ¬ 
ing upwelling pillars forming within bubbles. Our simulations show impacting 
Rayleigh-Taylor hngers reflecting on transition region plasma, ensuring that cool, 
dense chromospheric material gets mixed with prominence matter up to very 
large heights. This offers an explanation for the return mass cycle mystery for 
prominence material. Synthetic views at extreme ultraviolet wavelengths show 
remarkable agreement with observations, with clear indications of shear-flow in¬ 
duced fragmentations. 

Subject headings: magnetohydrodynamics (MHD) — Sun: filaments, prominences — 
Sun; corona 
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1. Rayleigh-Taylor in prominences 


Prominences are among ‘the most common features of the solar atmosphere’ (Parent! 


2014) and are a consequence of the Lorentz force levitating solar plasma against gravity. 


This creates density inversions in the hot solar corona favourable to Rayleigh-Taylor 
driven dynamics. Rayleigh-Taylor instability ~ the reason why water falls out of a cup 
turned upside down - can occur whenever a fluid or gas gets accelerated or pushed into 
a denser substance and is at the heart of many dynamical phenomena in astrophysical 
plasmas. Recent magnetohydrodynamic (MHD) modeling of Rayleigh-Taylor evolution for 
the Crab nebula showed a clear tendency to self-organize into larger-scale structures, with 


hlament sizes reaching up to a quarter of the entire pulsar wind nebula radius (Porth et al. 


2014a). Magnetic held-guided accretion processes onto magnetized stars are enriched by 


equatorially accreting Rayleigh-Taylor plasma tongues, protruding into the magnetosphere 


from the inner accretion disk edge (Kulkarni & Romanova 2008). In the solar context. 


Rayleigh-Taylor hlamentary structure can form during flux emergence and its reconnection 
with pre-existing coronal helds, as demonstrated by means of high resolution MHD 
simulations ( |Isobe et~ah 2005). 

Solar prominences also demonstrate Rayleigh-Taylor mediated dynamics, with 


Hinode Solar Optical Telescope observations (Berger et al. 2008) revealing how quiescent 


prominences show bright downflowing hlaments of several hundred kilometres in width, 
typical speeds of 0(10) km/s and ten minute lifetimes. At the same time, dark inclusions 
mark upflows at 20 km/s, rising up to 18 Mm heights and shedding voids that in turn grow 
to Mm sizes. Detailed observations showed how such dark upflows originate at the top 


of the chromosphere and can grow to 4-6 Mm plumes and rise to 15 Mm heights (Berger 


et al. 2010). They can form large-scale (20-50 Mm) buoyant arches or bubbles, and these 


rising bubbles were found to contain 25-120 times hotter material than the prominence 
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proper (Berger et al. 2011), strengthening the argnment for a magneto-thermal convection 
process typical for coronal cavity-prominence conhgnrations. Using a local model for 
a dipped prominence bottom bonndary, Rayleigh-Taylor mode development in three¬ 


dimensional (3D) MHD simulations demonstrated both upflows (Hillier et al. 2012a) and 
interchange reconnection leading to downward blob motions ( Hillier et al.||2012b ), in general 
agreement with observed local details. Recent modeling efforts have included partial 


ionization effects in local box models of the prominence-corona transition region (Khomenko 


et al. 2014), hnding clear differences with pure MHD approaches, as neutrals experience 


faster descents. Recent theoretical hndings quantihed the potentially stabilizing role of 
magnetic shear in idealized incompressible settings, important in the linear stages of 


Rayleigh-Taylor activity (Ruderman et al. 2014). A step towards global modeling of 
prominence dynamics in arcade systems conhrmed this role of sheared magnetic helds as 


well as the effects of line-tying on prominence stability (Terradas et al. 2015), but lacked 
the resolution to follow their development into the strongly nonlinear regime. In this paper, 
we set forth to realize this step, for the hrst time including chromosphere-transition region 
variations in a 3D prominence setup, and simulating well into the observed magneto-thermal 
convective motions. 


2. Numerical setup 

Our simulation box extends for 30 Mm horizontally (x) and vertically {y), and has a 
width (z) of 14 Mm. Using three levels of dynamic grid rehnement we achieve a resolution 
of 600 X 600 X 280, i.e. grid cell sizes down to 50 km. With the open-source MPI-AMRVAC 


software (Keppens et al. 2012 Forth et al. 2014b), we perform 3D, ideal MHD simulations 


for two cases that differ most markedly in the prevailing magnetic held strength throughout 
the domain: a weak held case at ~ 8 G, and a stronger held one at ~ 20 G. These 
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are representative values for quiescent prominence conditions. The initial magnetic field 
B = {Bx = 0.1998 G, 0, i? 2 (|/)) is purely planar and non-uniform, due to an exponential 
decrease of the strongest Bz{y) component in a layer of 2.5 Mm thickness above the initial 
prominence heigth i/p = 12.5 Mm. The analytic form for the initial Bz{y) is given by 

Bz = Bzo for y <yp, 

= Bzo exp [-{y - |/p)/Ab] for yp < y < yh 

= Bzo exp [-(i/b - |/p)/Ab] for yh < y ■ (1) 

The parameters are set to As = 15 Mm, and yh = 15 Mm, while Bzo differs between the 
weak and strong held case. This induces a local shear in magnetic held, and establishes 
an upward magnetic pressure force that lifts matter against solar gravity. The horizontal 
Bx held component leads to stabilization by tension forces against purely planar {x, y) 
Rayleigh-Taylor instabilities for all wavelengths exceeding about 33 km, slightly below 
our numerical resolution. Nonlinear ehects quickly dominate the dynamics at lengthscales 
fully captured in our study, a result corroborated by high resolution purely two-and-a-half 
dimensional simulations. The initial magneto-hydrostatic stratihcation introduces a 
transition region height at y^r = 2.5 Mm where the temperature smoothly connects an 
8000 K chromosphere to a 1.8 MK corona. 

Figure 1 shows in top and bottom left panels the temperature and density structure for 
the strongest magnetic held case. The dashed prohles above yp = 12.5 Mm quantify T{y) 
and p{y) exterior to the prominence, where the corona is isothermal but the density shows 
an increase due to magnetic levitation. Inside the prominence, the vertical stratihcation 
follows the solid curves shown in Figure 1: the prominence temperature is 12000 K in 
y G [yp, 15Mm], increases linearly with height to 60000 K at 25Mm, and there connects 
again to coronal temperatures above the prominence structure. This initial condition has 
the essential characteristics of solar hlaments, as the density contrast Pprom/Pcor below the 
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prominence is about 127.3 in this strong field case. When we set the overall dimensions 
of the prominence segment at 30 Mm length and 5 Mm width, we find that the total 
prominence mass is 7.5 x lO^^g for the weak field case, going up to 2.9 x lO^^g for the 
strong field case. These masses, together with the overall dimensions, all fall within their 
observationally known ranges. 


3. Global evolution 

The initial condition - though vertically in force balance ~ is out of pressure 
equilibrium in the x-direction across the prominence structure. This leads to a transient 
phase of successive compressions of the prominence matter (and in the coronal region 
above it), with shock waves traversing the periodic x-direction. These alter the detailed 
temperature-density variations throughout prominence and coronal regions upward from 
y = Upi but largely retain their essential characteristics, keeping the total prominence mass 
and typical corona-prominence density and temperature contrasts. More importantly for our 
study targeting long-term prominence internal dynamics, these transverse motions quickly 
become dominated by vertical and horizontal {x) velocity components, as demonstrated 
in the right panel of Figure 1, where the (scaled) kinetic energy evolution is plotted for 
each velocity component, for the strong field case (the weak field case behaves similarly 
in its energetic evolution). After about 4 minutes, vertical motions (solid line) dominate 
in kinetic energy, and they saturate before 10 minutes. Lateral movements (x-direction, 
dotted line) peak at about 11 minutes, while we ran our models for close to 15 minutes. The 
growth in vertical kinetic energy directly relates to Rayleigh-Taylor modes throughout the 
prominence segment, which are triggered by a superposition of 50 small-amplitude velocity 
perturbations that fit the periodic x-direction with random phases. Each individual flow 
perturbation represents a planar {vx,Vy,Q) incompressible velocity field, and is localized 
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about the bottom prominence position j/p and its midplane z = 0. 


Figure 2 gives a clear overview of the resulting prominence deformation and dynamics, 
by collecting a number of depth-integrated views taken at 6.9 minutes. This hgure is for 
the strong held case. Panels (a) and (c) provide views on the prominence when integrated 
along its length, showing its entire embedding within coronal plasma, while the nonlinear 
Rayleigh-Taylor development has created downward falling pillars that just reached 
transition region heights. Panel (a) integrates an additionally advected scalar, where green 
values correspond to prominence matter, dark purple indicates chromosphere plasma, and 
white is used for coronal material. Panel (c) relates to the instantaneous temperature 
structure, with white indicating cool (chromospheric and prominence) matter, and red 
is used for coronal values. This panel also shows a thin layer of hot matter just above 
the prominence structure, which locates shock-heated, initially evacuated matter found 
there. Animated views for the entire simulation in the format of Figure 2 are provided 
as online material, where the mentioned transverse compressions and their transient 
nature become evident. From our earlier simulations of actual prominence formation 
due to chromospheric evaporation, thermal instability and runaway catastrophic cooling 


events (Xia et al. 2011, 2012 Fang et ah 2013 Keppens & Xia 2014), these transient shock 


waves mimic the rebound shock fronts found to result from siphon flow driven impacts 
on the prominence-corona transition region. These rebound shocks ultimately repeatedly 
impact on the prominence structure, as a result of successive reflections when they reach 
chromosphere-corona transition regions along the heldlines. 


The edge-on views shown in panels (b) and (d) of Figure 2 contain the ^-integrated 
density structure, clearly dominated by falling Rayleigh-Taylor pillars with widths of up 
to 1000 km, and bubbles of upwardly curved prominence segments with lateral dimensions 
between 3000-4000 km. The resulting magnetic held deformation is visualized in panel 















(b), where streamlines, colored by the tracer from panel (a), are given for the 2 :-integrated 
in-plane magnetic held components. This shows how the falling pillars indeed interchange 
magnetic held structure, where we note that the prevailing plasma beta is typically 0.16 
(for the strong held case, and 0.38 for the weak one). The main displacements, as also seen 
in Figure 1, rapidly turn vertical and lateral, in accord with interchange modes that try to 
minimize held line bending, as the dominant magnetic held component is at all times. 
The same information can be deduced from panel (d) in Figure 2, where the superposed 
arrows likewise quantify the (; 2 -integrated in-plane) velocity structure. Clearly, regions 
with Rayleigh-Taylor hngers are overall downward-moving at this time, while hot coronal 
plasma shows the fastest upwelling hows within several of the bubbles. An interesting 
detail is the upwelling Rayleigh-Taylor hnger seen in panels (b)-(d) at horizontal distance 
X ~ 12.5 Mm that is seen to start at y ~ 10 Mm working its way upwards from the bottom 
region of a bubble. This bubble has just been closed from below, by the merging of two 
downwelling hngers that jointly continue their fall. A localized dense protrusion then 
swirls upwards, indicating that the relative acceleration (between light and dense matter) 
causing the Rayleigh-Taylor event now acts upwards in the bottom region of this bubble. 

A three-dimensional view on the prominence structure is given for about the same time 
in the left panel of Figure 3. This shows the temperature variation in vertical bounding 
planes at a; = 0 and = —7 Mm. It also shows an isosurface of the temperature marking 
the 30000 K isosurface, showing that all cool material is found in the downward pillars and 
at the lower regions of the bubbles. This isosurface also nicely traces out the location of the 
chromosphere-corona transition region, which has hardly been perturbed at this point in the 
evolution. The grey isosurface shows the rear-half of the tracer isosurface, at a value which 
locates the original prominence matter at all times, as well as the chromosphere-corona 
transition. In this view, we also see several of the upwelling Rayleigh-Taylor features, 
in the closed bubble discussed earlier but also near the x ~ 30 Mm front end. Hence 
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such temporary upwelling features with widths of about 500 km, should be identihable in 
the early stages of prominence formation and their internal dynamics. Note that speeds 
associated with individual larger-scale features, such as the falling and rising hlaments 
or bubbles, are several tens of km/s, also seen from the animated views provided. Using 
the tracer mentioned earlier, we quantify a prominence-material-only average vertical 
speed. This increases from zero up to about —20 km/s after 400 seconds, declining again 
afterwards. To quantify better actual speeds, we added to our 3D MHD simulation a set 
of 30 X 60 X 10 Lagrangian particles, which originally are positioned regularly on a grid 
throughout the simulation box. In the snapshot shown in Figure 3, twenty-four of the 
18000 particle trajectories obtained are visualized with streamlines that color from dark to 
white when time proceeds through the almost 15 minute interval simulated. These select 
initial locations all near the midplane x = 0, half of them initially right below the y = Up 
bottom prominence layer, while the other 12 are internal to the prominence. The right 
panel of Figure 3 shows the same set of particle trajectories, but viewed in {vy,y) phase 
space. Dotted horizontal lines mark the heights of the initial prominence edge yp, as well as 
the transition region height at 2.5 Mm. In this view, the ones that started internal to the 
prominence are indicated with larger symbols, while the thinner trajectories correspond to 
external (coronal) matter. Figure 3 shows that downward motions at up to 60 km/s prevail 
at hrst, but nearly all get deflected upwards after encountering the transition region. Both 
coronal and internal prominence matter can get accelerated up to speeds exceeding 120 
km/s. They can thereby reach heights well above their starting position, as some tracks 
go beyond 20 Mm height. Since a typical sound speed for the corona is 200 km/s, while 
the internal prominence sound speed is ten times lower, the process is highly nonlinear, 
and a vigorous magnetoconvective flow pattern extends from the transition region up into 
the prominence surroundings. Prominence matter can thereby repeatedly recycle, as it 
traverses a large range in altitude. These Lagrangian trajectories also imply that held lines 
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(mainly directed along z) indeed show significant interchange behavior. This aspect may 
be exaggerated in our simulation by the periodicity assumption: in reality, these field lines 
are part of an arcade system passing through the prominence matter, and line-tying effects 


play a role in determining their Rayleigh-Taylor stability properties (Terradas et al. 2015). 


An impression of the magnetoconvection that gets established throughout the 
prominence surroundings is shown in Figure 4, where the high field case is visualized at 
the endtime of our simulation, i.e. at 14.3 minutes. At left, we show the tracer isosurfaces 
that identify all prominence matter (colored by the local temperature), along with a 
grey isosurface that identifies the original chromosphere-corona transition region. The 
former isosurface shows that prominence bubbles have merged and grown into arch-shaped 
structures that can reach sizes up to 10 Mm in width. The latter isosurface demonstrates 
that the impacting Rayleigh-Taylor fingers can locally significantly perturb the transition 
region, and cause dense chromospheric matter to be hurled up to heights of 10 Mm or 
more. This provides an effective route to feed more cool, dense matter into the prominence 
environment, and hence plays an important role in its mass recycling. Figure 4 also shows 
the density structure in a vertical slicing plane at z = 0 in the box at right. In this view, 
we also visualize all tracer particles found between x = 15 and 30 Mm initially, where their 
color encodes the original starting height of the particle. At time zero, this color coding 
gives a plane-parallel green-orange-yellow-red distribution from top to bottom, while at 
the endtime from Figure 4, vigorous convection shows effective mixing in the entire region 
between 2.5 Mm and up to 23 Mm. While Figure 4 is for the high field strength case, we 
provide animated views for the low field strength case in the representation of Figure 4, 
as online material. Qualitatively similar trends occur in high and low field strength cases, 
although the maximal velocities attained are lower for the low field case, and the falling 
pillars reach the transition region a bit later in the evolution. This dynamical evolution 
allows one to interpret the temporal evolution of the component-wise kinetic energy through 
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the box shown in Fignre 1. The maximnm correlates with impacting falling pillars on the 
transition region, and lateral deflections maximize when np and down welling material meet 
np. 


4. Synthetic observations and conclusions 


Onr macroscopic simulations can be turned into extreme ultraviolet synthetic images, 


for direct comparison with those available from Solar Dynamics Observatory (Pesnell et al. 


2012) (SDO) observations using the Atmospheric Imaging Assembly (Lemen et ah 2012) 
(AIA). Especially its 304 A and 171 A channels provide views highlighting matter at 
80000 K and 800000 K, respectively. This then samples cooler prominence to transition 
region material. A synthetic view of both the high held (top rows) and low held case 
(bottom rows) in both EUV channels is given in Figure 5 at the endtime of our simulations, 
while animated views on the hnal seven minutes of evolution are provided online. One 
notices how cool prominence matter is found embedded in hotter material, with falling and 
rising structures over a fair range of lengthscales. The diherent wavelength channels show 
morphological diherences between hot and cool, up and downhow streams. The simulated, 
late nonlinear stages, especially for the low held case, show clear substructure developing 
along the edges of the largest bubbles, as seen in the bottom panels of Figure 5. At this 
stage, strong shear hows are established all along the arcs, that now extend as 10 Mm 
wide structures to heights of 18 Mm. We expect similar details to develop in the later 
stages of the high held strength case as well, as it also shows strong shear hows. This is in 
direct agreement with the latest observational details, pointing to Kelvin-Helmholtz and 


Rayleigh-Taylor interplay at the bubble boundary (Berger 2014). Visualizations of also 


the coronal channels (193 A and 211 A) further reveal the complex multi-temperature 
structure found in the magnetoconvective dynamics. Note that, by construction, our 
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side-on views show the prominence matter in emission, and assnme that the radiation is 
optically thin. This, together with the pure ideal MHD nature of our simulations, thereby 
neglecting important effects like coronal radiative losses, is an aspect to be improved upon. 
Further work can strive for even higher resolutions to capture smaller-scale hne-structure 
development from interplaying shear flow-driven, gravitational and thermal instabilities, or 
modihcations due to partial ionization conditions. Ultimately, ab initio simulations must be 


able to demonstrate the thermal instability mediated formation process of prominences (Xia 


et ah 2014), and simultaneously capture Rayleigh-Taylor mode development in realistic 


fluxrope-embedded prominence structures. 


This research was supported by the Interuniversity Attraction Poles Programme 
(initiated by the Belgian Science Policy Office, lAP P7/08 CHARM) and by the KU Leuven 
GOA/2015-014. Simulations used the VSC (Flemish Supercomputer Center) funded by the 
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Garching. C.X. is supported by FWO Pegasus hnancing. 
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Fig. 1.— The initial temperature (top left) and density (bottom left) stratification, both 
within (solid) and external (dotted) to the prominence. Right panel; the scaled kinetic 
energy evolution, plotted per velocity component: vertical (solid, y), lateral (dotted, x), and 
transverse (dashed, z). 
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Fig. 2.— Ray-traced views at about 6.9 minutes, along [panels (a), (c)] and across [panels 
(b), (d)] the prominence axis, showing contour views of: (a) the tracers used to identify 
prominence (green) and chromosphere (purple) material; (b)-(d) the density variation; (c) the 
temperature. Right top panel also shows heldlines based on integrated horizontal magnetic 
components, while bottom right arrows quantify the in-plane velocity variation. A movie in 
this format (for the high held case) is provided online. 
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Fig. 3.— Left: a 3D view on the prominence, at the same time as Figure 2, showing the 
temperature variation on bounding planes, as well as a (red) isocontour at 30000 K. The 
grey isosurface shows half of the prominence-bounding surface. Furthermore, 24 Lagrangian 
tracer paths are superposed, changing their color from black to white to indicate temporal 
variation. At right, the same 24 trajectories are displayed in a {vy,y) phase-space view. 
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Fig. 4.— After 14.3 minutes, this 3D view shows at left the temperature variation on the 
prominence boundary (in red to yellow), as well as (in grey) the location of the heavily 
perturbed chromosphere-corona transition region. The prominence is in a state of vigorous 
nonlinear magnetoconvection, also shown by its density variation in a cutting plane, and 
the tracer particles at right. The latter were originally arranged from green to red in plane- 
parallel fashion from top to bottom. While this hgure is for the high held case, a movie in 
this format for the low held case is provided as online material. 
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Fig. 5.— SDO AIA views on the endstate after 14.3 minutes for both the high held case 
(top row) and low held case (bottom row). Left panels are at 304 A, right panels for 171 A. 
A movie comparing both cases in this format is given online, covering the last 7 minutes of 


evolution. 



















